/*
 * Dmm_scf.cc
 *
 *  Created on: May 15, 2013
 *      Author: drewlewis
 */

#include "Dmm_scf.h"
#include "Density.h"
#include <TiledArray/algebra/diis.h>

using namespace sc;
namespace TA = TiledArray;

double DF_DMM(madness::World &world, ARRAY2 &D, ARRAY2 &S, ARRAY2 &H,
              ARRAY2 &F, ARRAY2 &G, ARRAY3 &Qpq, double nuc_repl){

    // create a DIIS object to extrapolate Fock matrix
    TA::DIIS<ARRAY2> diis;

    // How many purification steps to use in Density code
    int purification_steps = 2;
    int scf_iter = 1;
    double energy = 0;
    double error_norminf = 1.0;
    // Begin SCF iterations
    while(error_norminf > 4e-6){
        double iter0 = madness::wall_time(); // Iteration timer

        //Density Update
        double mad_conj0 = madness::wall_time();
        Density_Update(D, S, F, scf_iter, purification_steps);
        double mad_conjf = madness::wall_time();

        //Fock Build
        double Fock0 = madness::wall_time(); // Fock build timer
        G("i,j") =  2.0 * ( Qpq("i,j,X") * ( D("n,m") * Qpq("n,m,X") ) ) -
                        ( Qpq("i,n,X") * (  D("n,m") * Qpq("m,j,X") ) );
        F("i,j") = H("i,j") + G("i,j");
        double fence0 = madness::wall_time();
        world.gop.fence();
        double Fock_f= madness::wall_time(); // End fock build timer

        // Computing gradient for DIIS error calculation
        TA::Array<double,2> gradient = 8 * ( S("i,q") * D("q,x") * F("x,j") -
                                             F("i,q") * D("q,x") * S("x,j") );

        //Performing DIIS update of the Fock Matrix
        error_norminf = TA::expressions::norminf(gradient("i,j"));
        diis.extrapolate(F, gradient);

        //Energy for this iteration
        energy = TA::expressions::dot( 2.0 * H("i,j") + G("i,j"), D("i,j") );

        world.gop.fence(); // End of iteration work
        double iterf = madness::wall_time(); // End iteration timer

        //Printing a mess of stuff
        std::cout << "Iteration time = " << iterf - iter0 << "\n";
        std::cout << "Iteration(" << scf_iter <<"), energy = "
                        << std::setprecision(15) << energy + nuc_repl;
        std::cout << "  Error = " <<
                        error_norminf << std::endl;
        std::cout << std::setprecision(8);
        std::cout << "Fock build time = " <<
                        (Fock_f - Fock0) << std::endl;
        std::cout << "Fock fence time = " <<
                        (Fock_f - fence0) << std::endl;
        std::cout << "Density Upate time = " <<
                        mad_conjf - mad_conj0 << std::endl;
        std::cout << "\n\n";
        ++scf_iter;
    }

    //Final energy after convergence.
    std::cout << "Molecular energy = " << energy + nuc_repl << std::endl;

    return energy + nuc_repl;
}

double Direct_DMM(madness::World &world, ARRAY2 &D, ARRAY2 &S, ARRAY2 &H,
                  ARRAY2 &G, ARRAY2 &F,
                  TA::Array<double, 4, direct_tile_4c> &eri, double nuc_repl){
  // create a DIIS object to extrapolate Fock matrix
  TA::DIIS<ARRAY2> diis;

    // How many purification steps to use in Densiy code
    int purification_steps = 2;
    int scf_iter = 1;
    double energy = 0;
    double error_norminf = 1.0;
    // Begin SCF iterations
    while(error_norminf > 4e-6){
        double iter0 = madness::wall_time(); // Iteration timer

        //Density Update
        double mad_conj0 = madness::wall_time();
        Density_Update(D, S, F, scf_iter, purification_steps);
        double mad_conjf = madness::wall_time();

        //Fock Build
        double Fock0 = madness::wall_time(); // Fock build timer
        G("i,j") =  D("n,m") * ( 2 * eri("i,j,m,n")  - eri("i,n,m,j") );
        F("i,j") = H("i,j") + G("i,j");
        double fence0 = madness::wall_time();
        world.gop.fence();
        double Fock_f= madness::wall_time(); // End fock build timer

        // Computing gradient for DIIS error calculation
        TA::Array<double,2> gradient = 8 * ( S("i,q") * D("q,x") * F("x,j") -
                        F("i,q") * D("q,x") * S("x,j") );

        //Performing DIIS update of the Fock Matrix
        error_norminf = TA::expressions::norminf(gradient("i,j"));
        diis.extrapolate(F, gradient);

        //Energy for this iteration
        energy = TA::expressions::dot( 2.0 * H("i,j") + G("i,j"), D("i,j") );

        world.gop.fence(); // End of iteration work
        double iterf = madness::wall_time(); // End iteration timer

        //Printing a mess of stuff
        std::cout << "Iteration time = " << iterf - iter0 << "\n";
        std::cout << "Iteration(" << scf_iter <<"), energy = "
                        << std::setprecision(15) << energy + nuc_repl;
        std::cout << "  Error = " <<
                        error_norminf << std::endl;
        std::cout << std::setprecision(8);
        std::cout << "Fock build time = " <<
                        (Fock_f - Fock0) << std::endl;
        std::cout << "Fock fence time = " <<
                        (Fock_f - fence0) << std::endl;
        std::cout << "Density Upate time = " <<
                        mad_conjf - mad_conj0 << std::endl;
        std::cout << "\n\n";
        ++scf_iter;
    }

    //Final energy after convergence.
    std::cout << "Molecular energy = " << energy + nuc_repl << std::endl;

    return energy + nuc_repl;
}
